
#include "../include/Distributions.hpp"
#include "../include/utils.hpp"
//#include "../include/Gibbs.hpp"
//#include "../include/SmcSampler.hpp"
#include "../include/Resample.hpp"
#include "../include/HMC_trunc.hpp"
#include "../include/OrthantBS.hpp"
#include "../include/Orthan.hpp"
#include "../include/Orthan2.hpp"
#include "../include/Orthantst.hpp"
#include "../include/Orthanghk.hpp"
#include "../include/Orthantsampler.hpp"
#include "../include/Orthantsampler2.hpp"
#include "../include/Orthanghkcache.hpp"
#include "../include/Misc/EP_tg.hpp"
#include "../include/Misc/EP_tg_sep.hpp"
#include "../include/Orthanghkst.hpp"
#include "../include/Models/Thurston/Gibbs_Thurs.hpp"
#include "../include/Models/Thurston/PF_Thurstone.hpp"
#include "../include/Models/Thurston/PF.hpp"

extern "C"{
	void Gibbs(int *M, int *pX, int *p, int *n, double *res1, double *res2)
	{

		imat X=vect2mat<int,imat>(pX,*n,*p);
		Gibbs_Thurst GT(X);	
		vector<mat> B=GT(*M);
		mat2vect<double>(B[0],*M,*p,res1);
		mat2vect<double>(B[1],*M,1,res2);
	}
	void fPGibbs(int *M, int *pX, int *p, int *n, double *res1, double *res2)
	{

		imat X=vect2mat<int,imat>(pX,*n,*p);
		Resample::condMultinomial *R=new Resample::condMultinomial;
		PGibbs<Resample::condMultinomial> GT(X,R);	
		vector<mat> B=GT(*M,100);
		mat2vect<double>(B[0],*M,*p,res1);
		mat2vect<double>(B[1],*M,1,res2);
	}		
	
	void PF(int *M,int *Mx,int *pX, int *p, int *n, double *res1, double *res2)
	{

		imat X=vect2mat<int,imat>(pX,*n,*p);
		Gibbs_Thurst GT(X);	
		vector<mat> B1=GT(10);
		cout << X;
		mat m1=B1[0](99,span::all);
		double s1=B1[1](99,0);
		Resample::SysResample *R=new Resample::SysResample;
		PF_Thurstone<Resample::SysResample> PFT(X,R);	
		vector<mat> B=PFT(*M,*Mx,m1,s1);
		mat2vect<double>(B[0],*M,*p,res1);
		mat2vect<double>(B[1],*M,1,res2);
	}	

}
